################################################################
### This file contains code to do the synthetic control
### analysis in the paper "The Hirschman Effect of International
### Financial Ties", Xun Pang and Chong Chen, 2020
################################################################
setwd("Replication")

library("Synth")
load("SynthData.RData")

data2 <- SynthData
#tu <- unique(data2$statenme[data2$swap==1])

### Covariates for synthetic control analysis
controls <- c( "polity2_score", "sco_full", "tradedepc", "tradedepp", "pgdppc",
              "pgdpgrow", "dist", "partner_level", "bits_chn", "pta_china", "armtrade_log",
              "asean_memb", "atopally_us")

### Treated Units
treatno <- unique(data2$countryno[data2$swap==1])
allno <- 1:193

## Control units
controlno1 <- setdiff(allno, treatno)

## Get rid of units with missing data that disrupt the method
mis <- c()
for (i in controlno1){
    da <- subset(data2, countryno==i& year<2009)[,controls]
    ff <- sum(sapply(da, function(x)all(is.na(x))))
    kk <- sum(data2$distance_idealpoint_cn[data2$countryno==i])
    if(ff>0|is.na(kk)) mis <- c(mis, i)
    }
controlno <- setdiff(controlno1, mis) # 97 controls in total

Argno <- data2$countryno[data2$statenme=="Argentina"][1] # 34

tu <- unique(data2$statenme[data2$swap==1])

trno <- Argno

aa <-data2$swap[which(data2$countryno==trno)]
bb <- data2$year[which(data2$countryno==trno)]
adoptyr <- bb[aa==1][1] ## year of adoption

minyr <- min(data2$year)
maxyr <- max(data2$year)

## Data preparation
dataprep.out <-
              dataprep(foo = data2,
                       predictors = controls,
                       predictors.op = "mean" ,
                       time.predictors.prior = minyr:(adoptyr-1) ,
                       dependent = "distance_idealpoint_cn" ,
                       unit.variable = "countryno",
                       unit.names.variable = "statenme",
                       time.variable = "year",
                       treatment.identifier = trno,
                       controls.identifier = controlno,
                       time.optimize.ssr = (adoptyr-10):(adoptyr-1),
                       time.plot = (adoptyr-10):maxyr
                       )
	
### to see whether the variables should be rescaled
 dataprep.out$X1
 dataprep.out$Z1
	
### rescale some variables
dataprep.out$X1["pgdppc",] <- log(dataprep.out$X1["pgdppc",])
dataprep.out$X0["pgdppc",] <- log(dataprep.out$X0["pgdppc",])

dataprep.out$X1["dist",] <- log(dataprep.out$X1["dist",])
dataprep.out$X0["dist",] <- log(dataprep.out$X0["dist",])
	
### Synthetic control
 synth.out <- synth(data.prep.obj = dataprep.out,
                    method = "BFGS")
	
### See the gap
 gaps <- dataprep.out$Y1plot - (dataprep.out$Y0plot %*% synth.out$solution.w)
 gaps[1:3, 1]
	
### see the weights
 synth.tables <- synth.tab(dataprep.res = dataprep.out,
                           synth.res = synth.out
                           )
	
#synthTablesArg <- synth.tables
#save(synthTablesArg, file="synthTablesArg.RData")

couw <- synth.tables$tab.w
counw <- couw[order(couw$w.weights, decreasing = TRUE),]
## when do placebo test. Get 20 control units to do placebo test


### make plot  (Figure 5 in the artile)
 path.plot(synth.res = synth.out,
           dataprep.res = dataprep.out,
           Ylab = "Foreign Policy Distance to China",
           Xlab = "year",
           Ylim = c(0,2),
           Legend = c("Argentina","synthetic Argentina"),
           Legend.position = "topright"
           )
abline(v=2009,lty="dotted",lwd=2)

### make plot of the gap  (Figure 5 in the artile)
 gaps.plot(synth.res = synth.out,
           dataprep.res = dataprep.out,
           Ylab = "gap in real foreign policy distance with China",
           Xlab = "year",
           Ylim = c(-1.5,1.5),
           Main = NA
           )
abline(v=2008,lty="dotted",lwd=2)

argTE <- dataprep.out$Y1plot - (dataprep.out$Y0plot %*% synth.out$solution.w)
###################################################
### Placebo Tests
###################################################
## 1.  set 2005 to be the adoption year
## Control units
controlno1 <- setdiff(allno, treatno)
## Get rid of units with missing data that disrupt the method
mis <- c()
for (i in controlno1){
    da <- subset(data2, countryno==i& year<2005)[,controls]
    ff <- sum(sapply(da, function(x)all(is.na(x))))
    kk <- sum(data2$distance_idealpoint_cn[data2$countryno==i])
    if(ff>0|is.na(kk)) mis <- c(mis, i)
    }
controlno <- setdiff(controlno1, mis)# 97 controls in total

tu <- unique(data2$statenme[data2$swap==1])


## Argentina
trno <- data2$countryno[data2$statenme=="Argentina"][1] # 34
aa <-data2$swap[which(data2$countryno==trno)]
bb <- data2$year[which(data2$countryno==trno)]
adoptyr <- bb[aa==1][1]
minyr <- min(data2$year)
maxyr <- max(data2$year)
 dataprep.out <-
              dataprep(foo = data2,
                       predictors = controls,
                       predictors.op = "mean" ,
                       time.predictors.prior = minyr:2004 ,
                       dependent = "distance_idealpoint_cn" ,
                       unit.variable = "countryno",
                       unit.names.variable = "statenme",
                       time.variable = "year",
                       treatment.identifier = trno,
                       controls.identifier = controlno,
                       time.optimize.ssr = minyr:2004,
                       time.plot = minyr:maxyr
                       )

dataprep.out$X1["pgdppc",] <- log(dataprep.out$X1["pgdppc",])
dataprep.out$X0["pgdppc",] <- log(dataprep.out$X0["pgdppc",])

dataprep.out$X1["dist",] <- log(dataprep.out$X1["dist",])
dataprep.out$X0["dist",] <- log(dataprep.out$X0["dist",])

 synth.out <- synth(data.prep.obj = dataprep.out,
                    method = "BFGS")

 synth.tables <- synth.tab(dataprep.res = dataprep.out,
                           synth.res = synth.out
                           )

## when do placebo test. Get 20 control units to do placebo test

path.plot(synth.res = synth.out,
           dataprep.res = dataprep.out,
           Ylab = "Foreign Policy Distance to China",
           Xlab = "year",
           Ylim = c(0,2),
           Legend = c("Argentina","synthetic Argentina"),
           Legend.position = "bottomleft"
           )
abline(v=2004,lty="dotted",lwd=2)
abline(v=2008,lty="dotted",lwd=3)


### 2. Placeto Test on 20 control units

goodunit <-counw$unit.numbers[1:6]
store1 <- matrix(NA,length((adoptyr-10):maxyr),length(goodunit))
colnames(store1) <- goodunit

i <- 1
con <- goodunit
# run placebo test
for(iter in con)
 {
 dataprep.out <-
              dataprep(foo = data2,
                       predictors = controls,
                       predictors.op = "mean" ,
                       time.predictors.prior = minyr:(adoptyr-1) ,
                       dependent = "distance_idealpoint_cn" ,
                       unit.variable = "countryno",
                       unit.names.variable = "statenme",
                       time.variable = "year",
                       treatment.identifier = iter,
                       controls.identifier = controlno[controlno!=iter],
                       time.optimize.ssr = (adoptyr-10):(adoptyr-1),
                       time.plot = (adoptyr-10):maxyr
                       )

dataprep.out$X1["pgdppc",] <- log(dataprep.out$X1["pgdppc",])
dataprep.out$X0["pgdppc",] <- log(dataprep.out$X0["pgdppc",])

 dataprep.out$X1["dist",] <- log(dataprep.out$X1["dist",])
dataprep.out$X0["dist",] <- log(dataprep.out$X0["dist",])
	
# run synth
synth.out <- synth(
                   data.prep.obj = dataprep.out,
                   method = "BFGS"
                   )

# store gaps
 store1[,i] <- dataprep.out$Y1plot - (dataprep.out$Y0plot %*% synth.out$solution.w)
 i <- i+1
}

otherunit <- sample(counw$unit.numbers[-c(1:6)], 14)
store2 <- matrix(NA,length((adoptyr-10):maxyr),length(otherunit))
colnames(store2) <- otherunit

i <- 1
con <- otherunit
# run placebo test
for(iter in con)
 {
 dataprep.out <-
              dataprep(foo = data2,
                       predictors = controls,
                       predictors.op = "mean" ,
                       time.predictors.prior = minyr:(adoptyr-1) ,
                       dependent = "distance_idealpoint_cn" ,
                       unit.variable = "countryno",
                       unit.names.variable = "statenme",
                       time.variable = "year",
                       treatment.identifier = iter,
                       controls.identifier = controlno[controlno!=iter],
                       time.optimize.ssr = (adoptyr-10):(adoptyr-1),
                       time.plot = (adoptyr-10):maxyr
                       )

 dataprep.out$X1["pgdppc",] <- log(dataprep.out$X1["pgdppc",])
dataprep.out$X0["pgdppc",] <- log(dataprep.out$X0["pgdppc",])

 dataprep.out$X1["dist",] <- log(dataprep.out$X1["dist",])
dataprep.out$X0["dist",] <- log(dataprep.out$X0["dist",])
	
# run synth
synth.out <- synth(
                   data.prep.obj = dataprep.out,
                   method = "BFGS"
                   )

# store gaps
 store2[,i] <- dataprep.out$Y1plot - (dataprep.out$Y0plot %*% synth.out$solution.w)
 i <- i+1
}


# now do figure (Figure 5 in the article)
data <- data.frame(argTE, store1, store2)
rownames(data) <- (adoptyr-10):maxyr
colnames(data) <- c(goodunit, otherunit)

# Set bounds in gaps data
gap.start     <- 1
gap.end       <- nrow(data)
years         <- (adoptyr-10):maxyr
gap.end.pre  <- which(rownames(data)=="2008")

#  MSPE Pre-Treatment
#mse        <-             apply(data[ gap.start:gap.end.pre,]^2,2,mean)
#argentina.mse <- as.numeric(mse[1])
# Exclude states with 5 times higher MSPE than basque
#data <- data[,mse<100*argentina.mse]
Cex.set <- .75

# Plot
plot(years,data[,1],
     ylim=c(-1,1.5),xlab="year",
     xlim=c((adoptyr-10), maxyr),ylab="gap in foreign policy distance with China",
     type="l",lwd=2,col="black",
     xaxs="i",yaxs="i")

# Add lines for control states
for (i in 2:ncol(data)) { lines(years,data[gap.start:gap.end,i],col="gray") }

## Add Basque Line
lines(years,data[,1],lwd=2,col="black")

# Add grid
abline(v=2008,lty="dotted",lwd=2)
abline(h=0,lty="dashed",lwd=2)
legend("topleft",legend=c("Argentina","control countries"),
lty=c(1,1),col=c("black","gray"),lwd=c(2,1),cex=.8)
arrows(minyr,-1.5,(minyr+0.5),-1.5,col="black",length=.1)
text(2008.5,-1.5,"Signing BSA",cex=Cex.set)
abline(v=minyr)
abline(v=maxyr)
abline(h=-2)
abline(h=2)

### unit 35 ("Uruguay") and 101 ("Gabon") are lower than argentina
unique(data2$statenme[data2$countryno==35])


